Improved Parameterization of K + Production in p-Be Collisions at Low Energy 

Using Feynman Scaling 
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This paper describes an improved parameterization for proton-beryllium production of secondary 
K + mesons for experiments with primary proton beams from 8.89 to 24 GeV/c. The parameteri- 
zation is based on Feynman scaling in which the invariant cross section is described as a function 
of xp and pr- This method is theoretically motivated and provides a better description of the 
energy dependence of kaon production at low beam energies than other parameterizations such as 
the commonly used "Modified Sanford-Wang" model. This Feynman scaling parameterization has 
been used for the simulation of the neutrino flux from the Booster Neutrino Beam at Fermilab and 
has been shown to agree with the neutrino interaction data from the SciBooNE experiment. This 
parameterization will also be useful for future neutrino experiments with low primary beam energies, 
such as those planned for the Project X accelerator. 

PACS numbers: 13.25.Es,13.87.Ce 



I. INTRODUCTION 

This paper describes a parameterization for inclusive 
production of secondary K + mesons in proton-beryllium 
collisions, 



p + Be -> K+ + X. 



(1) 



for experiments with low primary proton beam energies 
ranging in kinetic energy from below 9 to 24 GeV. The 
parameterization is based on Feynman scaling(F-S) 
in which the invariant cross section is described as a func- 
tion of transverse momentum, px, and a scaling variable, 
xp = p„ /p„ M max . Various scaling parameterizations 
are known to describe data well above ~20 GeV 0, 
In this paper, we show that the F-S form describes data 
down to 8.89 GeV/c beam momentum. This result pro- 
vides an alternative model to the traditional "Modified 
Sanford-Wang" [H, Q parameterization used to describe 
secondary production at low primary proton beam mo- 
mentum. The results from this F-S analysis have been 
used in the neutrino flux parameterization of the Booster 
Neutrino Beam (BNB) at Fermilab and has been checked 
against measurements by the SciBooNE experiment @. 
This parameterization will be useful for future neutrino 
experiments using low primary proton beam energies. 

The primary motivation for this work was the simu- 
lation of neutrinos in the BNB line. This line provides 
neutrinos for the MiniBooNE Q and SciBooNE g exper- 
iments, as well as possible future experiments, including 
the upcoming MicroBooNE Q experiment. In this beam 
line, protons with 8 GeV kinetic energy are directed onto 
a 1.8 interaction length beryllium target. The charged pi- 
ons and kaons which are produced are focused by a mag- 
netic horn into a 50 m decay region, where they subse- 
quently decay to produce neutrinos. The average energy 
of 7r + (K + ) that decay to neutrinos in the MiniBooNE 
detector acceptance is 1.89 (2.66) GeV. Therefore 37.6% 



v„ from 7T & /i decoy 
y„ from K decay 
v, from 7T & /i decoy 
i/„ from K decay 




FIG. 1. Predicted and v a flux spectrum from decaying 
pions, kaons, and muons for the BNB and SciBooNE and 
MiniBooNE experiments. 



(92.1%) of the particles decay before the end of the 50 m 
long decay region. The most relevant decay modes for 
MiniBooNE are tt+ -> /i+i^, K + — > ii + v^, which pro- 
duce 99.4% of the neutrino beam, and K + — s- Tr°e + is e , 
fi + — > e + v^v e . l Kl — > 7r~ e + v e , and — >• ir + e~D e , which 
produce the remaining 0.6%. 

Figure [1] shows the predicted flux for the BNB line 
at the MiniBooNE detector. While the flux is predomi- 
nately due to 7r + decay, the K + decay is the dominant 
source above 2 GeV. The v e flux from kaon decay con- 
tributes one of the important backgrounds for neutrino 
oscillation searches looking for v e appearance. In addi- 
tion, the kaon neutrino flux provides an interesting source 
of high energy events for experiments on the BNB line 



for studying neutrino cross sections. Therefore, it is im- 
portant for the BNB line experiments to have a good 
first-principles prediction of K + production. 

A first-principles prediction for K + production is ob- 
tained from fitting data from secondary production ex- 
periments with primary beam momentum ranging from 
8.89 to 24 GeV/c. Nine data sets are considered, but 
only seven are used in the fit as it will be explained in 
Section IHII Because these data are taken at a range of 
beam energies, the data must be fit to a parameteriza- 
tion including changes with beam momentum in order 
to scale the result to the 8.89 GeV/c of the BNB line 
momentum. 
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Produced 


Exclusive 


M x 


\J ^thresh 


ijibeam. 
^thresh 


Hadron 


Reaction (GeV/c 2 ) 


(GeV) 


GeV 


7T + 


pn7T + 


1.878 


2.018 


1.233 


7T~ 


PP7T + 7T~ 


2.016 


2.156 


1.54 


7T° 


PP7T° 


1.876 


2.011 


1.218 


K+ 


A°pK+ 


2.053 


2.547 


2.52 


K~ 


ppK+K- 


2.37 


2.864 


3.434 


K° 


pE+K° 


2.13 


2.628 


2.743 



TABLE I. Threshold production channels for proton + pro- 
ton production of various mesons. The exclusive reaction is 
the final state with the minimum mass, Mx- y/sthresh and 
E'thr&sh are the threshold center of mass (CM) and laboratory 
energy. 



A. Feynman Scaling Formalism 

Over the past several decades, many experiments have 
made measurements of particle production by protons of 
various energies on many different nuclear targets. These 
data have been used to study the phenomenology of par- 
ticle production and have led to several scaling laws and 
quark counting rules. For inclusive particlcproduction, 
Feynman put forward a theoretical model |l| where the 
invariant cross section is only a function of xp and pp- 
The invariant cross section is related to the commonly 
used differential cross section by: 



d 2 a _ P^ E fv 
dpdn ~ E dp 3 ' 



(2) 



Defining 



d 3 a 



(3) 



this leads to: 



« G P at,, X 



(4) 



A is a factor and F is the F-S function that depends on 
xf and pp- The quantity p^ M max , which appears in the 
denominator of the definition of xp, depends upon the 
particle being produced and is derived from the exclusive 
channels given in Table fl] 

Feynman scaling has been demonstrated for secondary 
meson production at primary beam energies above ~15 
to 20 GeV this paper demonstrates the validity of 

F-S at lower primary beam energies for K + production. 
One might expect F-S to be a better parameterization of 
K + production than the "Modified Sanford-Wang" for- 
malism for two reasons. First, the F-S parameterization 
properly accounts for the kinematic effects of the large 
kaon mass where even at xp = 0, the outgoing kaon 
can have a significant laboratory momentum. Second, 
the functional form of the parameterization typically has 



peak production at xp — 0. This is in contrast to the 
"Modified Sanford-Wang" formalism, where the produc- 
tion rate continues to grow as xf becomes more negative. 



B. Feynman Scaling Parameterization for the 
Particle Production Cross Section 



The Feynman model can be used to describe the ex- 
pected xp and pr dependence using theoretically inspired 
functions for these dependences. For the xf depen- 
dence, a parameterization proportional to exp(— a I^fI^) 
or (1 — |£f|) c has the properties consistent with a flat 
rapidity plateau around xf = 0. The expectation of 
a limited pp range is provided by including exponential 
moderating factors for powers of pp- 

Using this guidance, a F-S parameterization has been 
developed to describe kaon production. In order to allow 
some coupling between the xp and pp distribution an 
additional exponential factor has been added that uses 
the product, \pp x xp\. The Cj's are the seven coefficients 
of the F-S function. The kinematic threshold constraint 

for K + production is imposed by setting — — — equal to 

dpdil 

zero for \xp\ > 1- 

Including these factors, the final parameterization has 
the form: 



dpdfl Ek 
x exp [c 3 \xf\ C4 



E 



K 



Pk 



d 3 <T 

dp 3 K ) \E K 
c 7 \p T x if| C6 



Ci x 



ClPT 



CbPr] 



(5) 



The "Modified Sanford-Wang" 
Parameterization 



Many neutrino experiments have used the "Modified 
Sanford-Wang" parameterization [H, Q (S-W): 
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dpdQ 



x exp 



ci PK C2 {l 



C 3PK 

Pbeam 05 



Pk 



PBEAM — Cg 

c§Qk(pk - c 7 pbeamcos C8 9k) 



(6) 



This functional form allows for some phenomenolog- 
ical parameterization of the variations associated with 
beam energy and process thresholds. As noted in one 
of the initial Sanford-Wang papers JH [H, the coeffi- 
cients for 7r + production are approximately given by: 
c-2 = 0.5, C4 = C5 = 1.67, and the cos# term is negligi- 
ble. With these substitutions, the formula shows a close 
although not perfect relationship with F-S (see Eq. [SJ 



E 



fa 
dp 3 



San ford— Wang 



A' F'(X) e 



-Cpi 



where 



and 



F'{X) = X l ' 2 {l-X)e 
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X 



V 



PBEAM 



(7) 



(8) 



(9) 



Therefore, the S-W fits to the K + data will show only 
approximate consistency with F-S. At low beam energy, 
produced particle mass effects can become important. 
Table U gives the minimum mass channels, their invari- 
ant mass, and the beam energy threshold for different 
particle production processes. In the S-W formula, the 
parameter eg is included to approximately provide the 
kinematic limit for produced particle momentum. Inves- 
tigations of the exact kinematic threshold for K + pro- 
duction show that the maximum pk is approximately 
equal to Pbeam — Pmff where Pmff varies from 1.7 
to 2.2 GeV as 9k goes from to 0.3 rad. One would 
therefore expect that eg would take on values similar to 



PDiff- On the other hand, the factor (1 



m \ 

A M — CO ' 



PBEAM— Cg ' 

introduces violations of the scaling behavior away from 
this limiting region. 

An additional problem with the S-W parameterization 
is that most of the function parameters (cj) will be ef- 
fectively fixed by the scaling constraints, and this will be 
limiting the flexibility of the function to match the xf 
and pt behavior. The parameter c 2 , for example, should 
be close to unity to provide the conversion from invariant 
to differential cross section. The parameter eg needs to be 
approximately equal to 2.0 GeV to provide the maximum 
Pk dependence, and the parameters C4 and C5 should be 
equal in order to preserve a basic xf dependence. Thus, 
the S-W parameterization has very little flexibility to fit 
the data distributions over the full kinematic range and 
therefore a formalism like Feynman scaling is required. In 
many of the following plots, we will compare prediction 
results coming from S-W and F-S parameterizations. 



II. EXTERNAL DATA SETS AND KINEMATIC 
COVERAGE 

Several K + production measurements have been made 
for beam momentum less than 25 GeV/c and are reported 
in Table HU Those experiments, except for Piroue, have 
beam momenta higher than the BNB value of 8.89 GeV/c 
although some of them such as Alcshin and Vorontsov are 
fairly close to the BNB beam momentum. The kaons that 
produce neutrinos in MiniBooNE span the kinematic re- 
gion with -0.1< x F <0.5 and 0.05< p T (GeV/c) <0.5 as 
shown in Figure [21 which is nicely covered by the experi- 
mental data sets listed in Table |TTJ Of course, we are us- 
ing the assumption that one can extrapolate these higher 
beam momentum data to the BNB energy value using a 
parameterization such as F-S. Thus, the first question to 
be answered is whether the data appears to follow these 
scaling parameterizations. 
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FIG. 2. Values of xf and pr for the data points of the various 
data sets in Table [TT1 The distribution for kaons that produce 
v e events in the MiniBooNE detector is shown as open boxes. 



The F-S hypothesis says that the invariant cross sec- 
tion should only depend on xf and pt- This hy- 
pothesis can further be tested by scaling all the data to a 
common beam momentum and checked by the behavior 
of the invariant cross section against the scaled value of 
Pk and 9k ■ Figure [3] shows the invariant cross section 
for scaled kaon momentum and angle bins using the F-S 
assumption. For this plot, the data from each data set is 
converted first to Xf and pt and then scaled to p K S9 and 
9 K 89 for a 8.89 GcV/c beam momentum. For example, 
given a cross section point at Pbeam=20 GeV/c with a 
given Pk and Ok , one can calculate the x f and pt for this 
point. One can then find the equivalent p' K and 9' K that 
would have the same Xf and pr at Pbeam=&-89 GeV/c. 
As seen from the plots, the data appears to obey the 
scaling hypothesis reasonably well except for the Lundy, 
Piroue, and Vorontsov data sets. Due to the disagree- 
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K+Data 


Ref. 


Pb (GeV/c) 


Pif(GeV/c) 


K (rad) 


xf 


p T (GeV/c) 


CTjVorm 


Abbott 


[101 


14.6 


2-8 


0.35-0.52 


-0.12 - 0.07 


0.2 - 0.7 


10% 


Aleshin 


[Hi 


9.5 


3 - 6.5 


0.06 


0.3 - 0.8 


0.2 - 0.4 


10% 


Allaby 


[121 


19.2 


3-16 


0-0.12 


0.3 - 0.9 


0.1 - 1.0 


15% 


Dekkers 


[13] 


18.8 , 23.1 


4-12 


0,0.09 


0.1-0.5 


0.0- 1.2 


20% 


Eichten 


[14] 


24.0 


4-18 


0-0.10 


0.1-0.8 


0.1 - 1.2 


20% 


Lundy 


[15] 


13.4 


3-6 


0.03,0.07,0.14 


0.1-0.6 


0.1 - 1.2 


20% 


Marmer 


[16] 


12.3 


0.5 - 1 


0,0.09,0.17 


-0.2- -0.05 


0.0-0.15 


20% 


Piroue 


[1JJ 


2.74 


0.5 - 1 


0.23,0.52 


-0.3- 1.0 


0.15-0.5 


20% 


Vorontsov 


[18] 


10.1 


1-4.5 


0.06 


0.03 - 0.5 


0.1-0.25 


25% 



TABLE II. Data sets for K + production with proton momentum lower 24 GeV/c. Pb indicates the beam momentum and o. 
gives the normalization error for the experimental data. 



merits of the Lundy and Piroue data, these data sets are 
not included in the fits described below. The Vorontsov 
data appears to agree in shape with the other data sets 
but has an anomalous normalization. Data sets not in- 
cluded in the fits are not discarded. They are compared 
separately to the fit results, as explained below. 



III. FEYNMAN SCALING AND 
SANFORD-WANG MODEL FITS TO THE K+ 
EXTERNAL DATA SETS 

Under the assumption that the experimental data fol- 
low the Feynman or S-W scaling models, we can deter- 
mine a parameterization that best fits these data sets. 
The various production data sets are used as input to a 
fit for the scaling function parameters that best describe 
the data. The fit uses a x 2 minimization technique using 
Minuit [III to perform the numerical minimization. Each 
experiment is allowed to have an independent normaliza- 
tion parameter that is constrained by the published nor- 
malization uncertainty. The fit minimizes the following 
function for an experiment j: 



given in Eq. [5] (or |6l and for the normalization factors 



x 2 = 



E 



(Nj x SF t - Data t ) 2 
(f x °if 



(l-#i) S 



,(10) 



'JV, 



where i is the (Pk$k) bin index, SF is the scaling func- 
tion prediction evaluated at the given (j>beam,&k ,Pk), 
Datdi is the measurement at a given (pbeam,@k,Pk), 
Oi is the data error for measurement i, f is the scaling 
factor to bring the \ 2 /d.o.f. = 1, Nj is the normalization 
factor for experiment j, is the normalization uncer- 
tainty for experiment j, and d.o.f. indicates degree of 
freedom. The total \ 2 for external data sets is then the 
sum over the experiments of the individual xj values, 



(ii) 



The x 2 is minimized in order to obtain the best values 
and uncertainties for the parameterization coefficients Cj , 



N, 



The uncertainties on the fit values at la are de- 



2 

min 



termined from a Ax 2 = 1 change with respect to x 
and the fit also yields a covariance matrix that can be 
used to propagate correlated errors associated with the 
parameterization of the cross section. 

A F-S fit to all the experimental data sets with 0.0< 
P£ 89 (GeV/c) <6.0 gives a x 2 /d.o.f. equal to 4.03 with 
large x 2 contributions from data with P|-- 89 <1.2 GeV/c 
and P|P 9 >5.5 GeV/c. Therefore for the final scaling 
fits, the points with the larger pull terms, defined as 
((Nj x SFi — Datdi) /(Ti), have been eliminated by only 
using data with 1.2< Pf < 89 (GeV/c) <5.5. 

The 1.2 GeV/c cut effectively removes data at nega- 
tive Xf where the nuclear environment starts to play an 
important role. This cut also eliminates all the Marmer 
data points. 

With all of these requirements, the x 2 /d.o.f. for the 
F-S fit is reduced to 2.28. The uncertainties for the fitted 
cross section need to be corrected for this x 2 /d.o.f., which 
is larger than 1.0. This is accomplished by scaling up the 
errors of each of the data points by sJx 2 /d.o.f . before 
doing the fit. Figure |4] shows the pull terms for the seven 
parameter F-S fit where the errors have been scaled up 
by this v / X 2 /d.o.f. = 72^28 =1.51. 

A S-W fit to all the experimental data has been per- 
formed as well. To be able to directly compare the S-W 
with the F-S fit we have included in the S-W fit only 
data with 1.2< P| 89 (GeV/c) <5.5. The x 2 /d.o.f. for 
the S-W fit is equal to 6.05. 



IV. COMPARISON OF FEYNMAN SCALING 
TO SANFORD WANG RESULTS AND 
NEUTRINO PREDICTIONS 



Tables IXTll and ITVl report the final fit values for the co- 
efficients and the normalization factors for the F-S and S- 
W parameterizations respectively. Figures [5] and |5] show 
the fit function curves for the F-S and S-W parameter- 
izations as compared to the data. The fits are stable 
with respect to parameter starting values and yield pos- 
itive definite covariance matrices. The error bands in 
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FIG. 3. K + production data sets scaled to the MiniBooNE beam momentum of 8.89 GeV/c using F-S. The Y-axis units are 
(mb x c 3 /GeV 2 ). The production angle varies from to 0.225 rad. 



Figures [5] and [6] are determined by propagating the co- 
variance matrix for the Cj parameters to the invariant 
cross section errors. 

As seen from the plots, the F-S function gives a very 
good description of the data over the full kaon momentum 
range used in the fit and has a reasonable x 2 /d.o.f. = 2.3. 
Below 1.2 GcV/c, the F-S prediction has some disagree- 
ment with a few the of the Marmer (not included in the 
fit) and Abbott data points but in general is also fitting 
that well in that region. The normalization factors for the 
F-S fits are within la of the quoted experimental error 
except for the Vorontsov data (see Table HIT)) . As men- 
tioned above, the Vorontsov data shows a systematically 
low normalization with respect the other sets of scaled 
data. Therefore, for all the scaling fits, the Vorontsov 
data has only been used for shape information by giving 
the normalization a large uncertainty (500%). 

In contrast, the S-W final fit parameterization has 



rather large discrepancies with the data in almost all re- 
gions and has a much larger \ 2 /d.o.i. = 6.05. Addition- 
ally, the normalization factors given in Table |IV] are very 
much outside of the quoted experimental errors and, for 
example, the factors for Eichten and Allaby differ from 
1.0 by 2 to 3cr. 

Tables fVl and IVT1 list the differential cross sections for 
several different kinematic points for kaon production. 
The uncertainties are obtained by propagating the co- 
variance matrix for the Cj coefficients into the scaling 
function. The first three points in Table [V] and IVII cor- 
respond to the mean kaon production points that pro- 
duce electron neutrino of 0.35, 0.65, and 0.95 GeV in 
MiniBooNE. The fourth point corresponds to the kaon 
kinematics that produce average energy neutrinos from 
all kaon decays (called the "kaon sweet spot"), and the 
fifth point is associated with the mean kaon kinematics 
for the highest energy kaon-decay muon neutrinos ob- 



G 




Sigma 



FIG. 4. Values of the pull terms, (Nj x SFi - 
Datai) / Oi, for each data point for the F-S fit for 1.2< 
Pfc S9 (GeV/c) <5.5. The data errors, <7i, have been scaled 
up by \/x 2 /d.o.f '. — V2.28 = 1.51. The Gaussian fit gives a 
X 2 /number of degree of freedom — 35.51/35, with a mean 
value = (-0.18 ± 0.11) and sigma = (0.90 ± 0.17). 



Feynman Scaling 


1.2< P K sa (GeV/c) <5.5 




Fit 


Value 


Error 




cl 


11.70 


1.05 




c2 


0.88 


0.13 




c3 


4.77 


0.09 




c4 


1.51 


0.06 




c5 


2.21 


0.12 




c6 


2.17 


0.43 




c7 


1.51 


0.40 


Input Error 


Aleshin 


1.09 


0.07 


0.10 


Allaby 


1.04 


0.07 


0.15 


Dekkers 


0.84 


0.06 


0.20 


Vorontsov 


0.53 


0.04 


5.00 


Abbott 


0.76 


0.07 


0.15 


Eichten 


1.00 


0.07 


0.15 


X 7d.o.f. (no /) 


2.28 


(d.o.f. = 119) 





TABLE III. Results for the F-S fits to the K + data including 
a single normalization factor for each experiment. The data 
errors have been scaled up by a factor of \J x 2 1 'd.o.f '. ={= 1.51 
when included in the fit but the x 2 /d.o.f. value listed is for 
the data without this scaling, d.o.f. indicates here degree of 
freedom and "no /" means no correction factor applied. 



served in MiniBooNE. As seen from Table [V] and I VII , 
the two parameterizations give much different results for 
the cross section values and uncertainties with the F-S fit 
giving a larger value by a factor 2 for the lowest energy 
neutrino bin at 0.35 GeV. The source of the this discrep- 
ancy is a large drop in the invariant cross section of the 
S-W parameterization at large angles. 



Modified S-W 


1.2< P K 


89 (GeV/c) <5.5 




Fit 


Value 


Error 




cl 


14.89 


1.89 




c2 


0.91 


0.13 




c3 


12.80 


7.46 




c4 


2.08 


0.35 




c5 


2.65 


0.50 




c6 


4.61 


0.10 




c7 


0.26 


0.01 




c8 


10.63 


7.06 




c9 


2.04 


0.01 


Input Error 


Aleshin 


1.02 


0.09 


0.10 


Allaby 


0.74 


0.09 


0.15 


Dekkers 


0.57 


0.08 


0.20 


Vorontsov 


0.42 


0.04 


5.00 


Abbott 


1.38 


0.11 


0.15 


Eichten 


0.59 


0.08 


0.15 


X7df (no /) 


6.05 


(d.o.f. = 117) 





TABLE IV. Results for the S-W scaling fits to the K + 
data including a single normalization factor for each exper- 
iment. The data errors have been scaled up by a factor 
of -y/x 2 /d.o.f '. =f= 2.46 when included in the fit but the 
X 2 /d.o.f. value listed is for the data without this scaling, 
d.o.f. indicates here degree of freedom and "no /" means no 
correction factor applied. 





r>8.89 

(GeV/c) 


0k 
(rad) 


&K prod 

(mb) 


E v = 0.35 GeV 


1.52 


0.213 


9.37 ± 0.73 (7.8%) 


E v = 0.65 GeV 


2.07 


0.127 


10.69 ± 0.75 (7.0%) 


E v = 0.90 GeV 


2.45 


0.103 


10.22 ± 0.71 (6.9%) 


Kaon Sweet Spot 


2.80 


0.106 


8.67 ± 0.60 (6.9%) 


HE Events 


4.30 


0.055 


4.73 ± 0.33 (7.0%) 



TABLE V. Differential cross section values for various kine- 
matic points for the 1.2< P K <5.5 GeV/c F-S fit. The first 
three results are for the average kaon kinematics that give 
electron neutrinos with the given energy. The last result is 
for the average kaon kinematics associated with highest en- 
ergy j/ m events in MiniBooNE. 





D 8.89 


0k 


0~K prod 




(GeV/c) 


(rad) 


(mb) 


E v = 0.35 GeV 


1.52 


0.213 


4.25 ± 0.77 (18%) 


E v = 0.65 GeV 


2.07 


0.127 


8.99 ± 1.34 (15%) 


E v = 0.90 GeV 


2.45 


0.103 


9.91 ± 1.43 (14%) 


Kaon Sweet Spot 


2.80 


0.106 


7.73 ± 1.13 (15%) 


HE i/ M Events 


4.30 


0.055 


5.24 ± 0.84 (16%) 



TABLE VI. Differential cross section values for various kine- 
matic points for the 1.2< Pk <5.5 GeV/c S-W scaling fit. 
The first three results are for the average kaon kinematics 
that give electron neutrinos with the given energy. The fourth 
result is the previous point used for a kaon sweet spot. The 
last result is for the average kaon kinematics associated with 
highest energy events in MiniBooNE. 
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FIG. 5. Invariant kaon production cross section in mb x c 3 /GeV 2 versus kaon momentum for all data along with the results 
of the F-S fit to data with 1.2< P^ 89 GeV/c <5.5. The Pk, Ok, and invariant cross section fits and the data points have been 
scaled to a beam momentum of 8.89 GeV/c assuming F-S and normalized according to the fit results. This plot shows data and 
fit results for various value of 8 in bins from to 0.225 rad. The three solid curves show the central value and la uncertainty 
for the F-S fit. 



The predictions for the size and kinematic dependence 
of the invariant differential cross section as function of 
K + momentum arc quite different for the F-S and S-W 
parameterizations as shown in Figure[71 especially for low 
value of the K + momentum. 



To illustrate the difference between the F-S and the 
S-W predictions, we have used an analytic simulation 
of the BNB neutrino beam line designed for the Mini- 
BooNE experiment (described in Reference |2(j). Ta- 
ble IVIII gives the comparison of the predicted v e event 
rate from K + — > i:e + v e using the above F-S and S-W pro- 
duction parameterizations as calculated using this BNB 
simulation. 



V. HIGH ENERGY PARAMETERIZATION 

The hypothesis of F-S has also been verified to hold 
with different parameterizations over a wide range of pri- 
mary proton beam energies (from 24 GeV to 450 GeV) . In 
Bonesini et al. data at higher proton energies has been 
empirically parameterized as a function of the transverse 
momentum (j>t) and the scaling variable xr = E* / E* nax 
where E* is the energy of the particle in center-of-mass 
frame. The choice of these variables for the description 
of the invariant cross section (radial scaling) is motivated 
again by an assumed scaling behavior of the invariant 
cross section. The radial scaling variable is approxi- 
mately equal to the F-S variable at high energy and has 
the property of never taking on a negative value. (A dc- 
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P* (GeV/c) 



FIG. 6. Invariant kaon production cross section in mb x c 3 /GeV 2 versus kaon momentum for all data along with the results of 
the S-W scaling fit to data with 1.2< P^ 89 GeV/c <5.5. The Pk,9k, and invariant cross section fits and the data points have 
been scaled to a beam momentum of 8.89 GeV/c assuming F-S and normalized according to the fit results. This plot shows 
data and fit results for various value of 9 in bins from to 0.225 rad. The three solid curves show the central value and ler 
uncertainty for the S-W scaling fits. 



9 K 




Feynman Scaling Fit 




K 3 


Sanford-Wang Fit 




Angular Bins (rad) 


All E„(GeV) 


<1 GeV 


>2 GeV 


All £„(GeV) 


<1 GeV 


>2 GeV 


0.015 


36.7 


2.6 


18.0 


43.4 


3.3 


19.3 


0.045 


92.5 


8.4 


35.9 


111.0 


12.0 


35.9 


0.075 


110.5 


13.7 


27.0 


141.3 


22.6 


26.5 


0.105 


96.8 


17.2 


4.4 


138.3 


32.6 


4.1 


0.135 


59.1 


21.8 


0.0 


100.5 


45.8 


0.0 


0.175 


39.4 


32.4 


0.0 


83.7 


73.9 


0.0 


0.225 


21.9 


21.9 


0.0 


56.8 


56.8 


0.0 


Total 


476.6 


137.9 


85.3 


731.2 


303.4 


85.9 



TABLE VII. Electron neutrino event rate in MiniBooNE for 5.0x10 proton on target for K+ 3 decays with F-S and S-W 
parameterizations. The events were calculated using MiniBooNE simulation and are for a beam radius less than 6.0 m. The 
different columns list the selected electron neutrino events for all E v , E v < 1 GeV, and E v > 2 GeV. Uncertainty in the 
neutrino event rate due to the F-S or S-W parametrization is 7% and 15% respectively as described in Table IVl and I VII 
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FIG. 7. Invariant kaon production cross section in units of mb x c 3 /GeV 2 versus kaon momentum in GeV/c for the S-W, F-S 
and radial scaling (Bonesini)Q] parameterizations for a beam momentum of 8.89 GeV/c. The results are shown for various 6 
bins from to 0.225 rad. The three solid curves, respectively for the F-S and S-W fits, show the central value and la uncertainty 
for each of the fits. 
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tailed comparison of radial scaling and F-S can be found 
in 0, [2l| , where the authors compare different models 
with the production data at different energies down to 
about 24 GeV.) 

Boncsini et a/.0] has obtained an empirical parame- 
terization based on radial scaling fits to data collected 
with 400 GeV/c and 450 GeV/c protons incident on a 
Be target. The results from this parameterization are 
compared in Figure [7] to the predictions of F-S and S-W 
models at a proton momentum of 8.89 GeV/c. As seen 
from Figure [71 this radial scaling model underestimates 
K + production at a beam momentum of 8.89 GeV/c by 
more than a factor of two even though the parameter- 
ization describes well the high proton momentum data 
(>24 GeV/c) 0. 

VI. THE SCIBOONE K + MEASUREMENTS 

The SciBooNE collaboration has reported a measure- 
ment [22[ for K + production in the BNB with respect to 
the Monte Carlo (MC) beam simulation. The SciBooNE 
experiment collected data in 2007 and 2008 with neutrino 
(0.99 x 10 20 protons on target (POT)) and antineutrino 
(1.53 x 10 20 POT) beams in the FNAL BNB line. The 
SciBooNE detector is located 100 m downstream from 
the neutrino production target. The flux- averaged mean 
neutrino energy is 0.7 GeV in neutrino running mode and 
0.6 GeV in antineutrino running mode 

The SciBooNE detector consists of three detector com- 
ponents; SciBar, Electromagnetic Calorimeter (EC) and 
Muon Range Detector (MRD). SciBar is a fully active 
and fine grained scintillator detector that consists of 
14,336 bars arranged in vertical and horizontal planes. 
SciBar is capable of detecting all charged particles and 
performing dE/dx-based particle identification. The 
EC is located downstream of SciBar. The detector is 
a spaghetti calorimeter with thickness of 11 radiation 
lengths and is used to measure ir° and the intrinsic u e 
component of the neutrino beam. The MRD is located 
downstream of the EC in order to measure the momen- 
tum of muons up to 1.2 GeV/c with range. It consists 
of 2-inch thick iron plates sandwiched between layers of 
plastic scintillator planes. 

In the SciBooNE experiment, particle production is 
simulated using the methods described in Ref. [2(| ■ The 
production of K + is simulated using the F-S formalism 
as described in Section iTAl with the coefficients reported 
in Table IIIII The predicted double differential cross sec- 
tion at the mean momentum and angle for kaons which 
produce neutrinos in SciBooNE (px = 3.87 GeV/c and 
9 K = 0.06 rad) is 



p K and K is 7%. The SciBooNE and MiniBooNE collab- 
oration have adopted a conservative error of 40%. This 
larger error was chosen because of the uncertainties in 
extrapolating the K + prediction data from high to low 
proton beam energy using the F-S and S-W models as 
explained in References [20, [HJ . 



A. SciBooNE K + Production Measurement 

The SciBooNE data can be used as an additional con- 
straint in fits to K + production cross sections. In Sci- 
BooNE, neutrinos from K + decay are selected using high 
energy interactions within the volume of the SciBar 
detector. The high-energy selection is accomplished by 
isolating charged current interactions that produce a 
muon that crosses the entire MRD. This sample is further 
divided into three sub-samples based on whether 1, 2, or 
3 reconstructed SciBar tracks are identified at the neu- 
trino interaction vertex in the SciBar detector. Since the 
reconstruction of the energy of the muon is not possible 
because the muon exits the MRD detector, the recon- 
structed muon angle relative to beam axis is used as the 
primary kinematic variable to separate neutrinos from 

pion and kaon decay. The values for — — — - for neutrino, 

dpdil 

antineutrino, and combined data mode running are given 
in Table I Villi along with the mean energy and angles for 
the corresponding K + samples. The F-S and S-W pre- 
diction values are obtained using the parametrizations 
described in Section II Al and II CI along with the parame- 
ters listed in Table HTT1 and IP71 





E K+ (GeV) 


6 K+ (rad) 


(fa 

dpdQ 
(mb/(GeV/c x sr) 


^-mode 


3.81±0.03 


0.07±0.01 


5.77±0.83 


u-mo&e 


4.29±0.06 


0.03±0.01 


3.18±1.94 


v + i/-mode 


3.90±0.03 


0.06±0.01 


5.34±0.76 


F-S prediction 


3.90 


0.06 


6.30±0.44(7%) 


S-W prediction 


3.90 


0.06 


6.84±1.09(16%) 



TABLE VIII. Measured — — — , mean energy, and mean angle 
dpdil 

(with respect to proton beam direction) for the selected K + 
in neutrino, antineutrino, and the combined neutrino and an- 
tineutrino samples using MiniBooNE MC. Errors on the mean 
energy and mean angle values correspond to the error on the 
mean for the relative distributions. F-S and S-W predictions 
are also reported at the mean SciBooNE K + energy and an- 
gle. 



— — = (6.3 ± 0.44) mb/{GeV/cx sr), (12) 
dpdil 

The error on the double differential cross section pre- 
diction using the F-S parametrization at the SciBooNE 



The K + momentum versus angle distribution for the 2- 
track SciBar sample in the simulation is shown in Fig- 
ure m 

Figure [5] shows the kinematics of the selected K + 
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□> 
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■o 
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0- 
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K + Production Momentum (GeV/c) 



1 



(a) Neutrino Sample 




K + Production Momentum (GeV/c) 



(b)Antineutrino Sample 

FIG. 8. The true K + momentum versus angle distribution in 
the SciBooNE MC for neutrino-mode (on top) and antineu- 
trino mode (on bottom) running. The unit for the color scale 
is number of events POT normalized. 



The SciBooNE measurement is a direct test of the ex- 
trapolation of parameterizations found from higher beam 
energies to the MiniBooNE beam energy. The predic- 
tions for the double differential cross section for the F-S 
and S- W models are reported in Table I Villi and shows a 
good agreement with the SciBooNE measurement, a bet- 
ter agreement is found in the case of the F-S paramctriza- 
tion. 

The SciBooNE K + production measurements can also 
be added to the F-S fit as additional external data using 
the following procedure. First we retrieve all the Sci- 
BooNE MC K + events with their Oi and pi for the neu- 
trino and antineutrino sample. Then we calculate the 
following quantities: 



d 2 G 

dpd£l 



{cf it ,0i,pi) 



d 2 a 
dpdQ 



{cAIC,0i,Pi) 



(13) 
(14) 



These quantities are then used at each fit step to build 
a pull term, defined in Eq. [151 to be added to the % 2 of 
the fit. 



pull — term v ^ 



(Hl 



Kprod,SB 



error K+ od SB 



(15) 



events in SciBooNE, while Figure H] shows the kinemat- 
ical region as function of angle and momentum for K + 
mesons that produce v e events in MiniBooNE. 




P„ (GeV/c) 



FIG. 9. Kinematical region as function of angle and momen- 
tum for the K + mesons that produce v e events in MiniBooNE. 
The unit for the color scale is number of events. 



Each data point in 9j and Pi is rewcightcd using the 
double differential cross section value for the current set 
of 

Ci coefficient of Eq. [5] computed at each step of the Mi- 
nuit fit. The set of coefficient used in the MC is labeled 
as cmc, the values of these coefficients is listed in Ta- 
bleHIIl The K+ od SB and error K+ B in the Eq.Hg] 
are the values of the SciBooNE production measurement 
and error (see Table lyTn)) . respectively. 

Two separate pull terms are added to the fit % 2 cor- 
responding to the SciBooNE neutrino and antineutrino 
K + production measurements. 

The results of scaling function fit to all experiments 
with 1.2< PjP 9 <5.5 GeV/c, including the SciBooNE 
data, are given in Table IIXI The covariance matrix is 
given in Table [X] and the correlation matrix is presented 
in Figure [TOl 

Table KJ lists the differential cross section for the kaon 
production at the various kaon kinematic points. The 
uncertainties are obtained as described in Section [IVI 

Tabic |X] gives the covariance matrix for the base- 
line scaling fit using kaon production data with 1.2< 
p8;89 <5_5 GeV/c. The correlation matrix is basically 
made of two blocks, one associated with the c\ through 
C7 parameters and one associated with the experimental 
normalization factors. The only coupling of these two 
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Scaling Fits 


1.2< Pi- 
Value 


S9 (GeV/c) <5.5 
Error 




cl 


11.29 


0.93 




c2 


0.87 


0.13 




c3 


4.75 


0.09 




c4 


1.51 


0.06 




c5 


2.21 


0.12 




c6 


2.17 


0.43 




c7 


1.51 


0.40 


Input Error 


Aleshin 


1.12 


0.07 


0.10 


Allaby 


1.07 


0.06 


0.15 


Dekkers 


0.87 


0.06 


0.20 


Vorontsov 


0.55 


0.04 


5.00 


Abbott 


0.79 


0.07 


0.15 


Eichten 


1.03 


0.06 


0.15 


X7df (no /) 


2.28 


(d.o.f. = 119) 





TABLE IX. Results for the F-S fits to the K + data including 
a single normalization factor for each experiment and includ- 
ing the two SciBooNE pull term constraints. Error treatment 
is the same as described in Section ITiTl d.o.f. indicates here 
degree of freedom and "no /" means no correction factor ap- 
plied. 



sets is through c\ which has significant correlations with 
the normalization factors. This is expected since the c\ 
parameter sets the normalization of the scaling function 
and should be determined by the data normalizations. 



>.0.125 

I 0.12 

g0.115 

3 0.11E- 
u 

£0.105 

CO 

I 0.1 
0.095 

0.09 E- 
0.085, 



-Without SB in Fit 
-With SB in Fit 



0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 

Angle (rad) 



FIG. 11. Relative uncertainty on the double differential cross 
section as function of K + angle (0.0< 0k <0.25 rad) pre- 
dicted by the F-S with and without including the SciBooNE 
production measurement. 



sJ>.125 
'J 0.12 
g0.115 
« 0.11 
!o.105 

DC 

0.1 
0.095 
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FIG. 12. Relative uncertainty on the double differential cross 
section as function of K + momentum (1.2< Pk <5.5 GeV/c) 
predicted by the F-S with and without including the Sci- 
BooNE production measurement. 




FIG. 10. Correlation for the seven parameters in the F-S fit 
function and six normalization factor parameters after apply- 
ing the SciBooNE constraint to the fit due to the K + produc- 
tion measurement. 

The terms of the covariance matrix from the F-S fit 
that includes the SciBooNE production measurement in- 
clude the factor 1.51 for the data set errors rescaling. 

The relative uncertainties on the predicted double dif- 
ferential cross section by the F-S fit as function of K + 
angle and momentum decrease including the SciBooNE 
measurement as shown in Figures [TT] and 1121 



The SciBooNE measurement confirms the validity of 
the F-S parametrization and including the SciBooNe 
measurement as an additional experimental data to the 
Feynman Scaling fit contributes in improving both the 
error uncertainty on the parametrization coefficients and 
in lowering the total uncertainty in the predicted K + 
production at 8.89 GeV/c proton momentum. 



B. SciBooNE K + Rate Measurement 

In addition to a measurement of K + production, the 
SciBooNE collaboration has also published a measure- 
ment of the observed to MC predicted ratio for K + 
produced neutrinos and antineutrinos interacting in the 
SciBar detector. The results are summarized in Ta- 
ble EU The SciBooNE rate is the product of the K+ 
production and neutrino cross section on carbon as ex- 
plained in Ref. [22[. Since this result also includes the 
neutrino interaction cross section, it cannot be directly 
compared with the other experimental data presented in 
Table|n] This constraint not only covers the neutrino flux 
from K + decay but also constrains the neutrino interac- 
tion cross section because the two targets are composed 
of similar material. It should be noted that this analy- 
sis is a specific application to MiniBooNE and is not a 
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cl 


c2 


c3 


c4 


c5 


c6 


c7 


cl 


0.84 


0.48E-01 


0.39E-02 


-0.32E-01 


-0.36E-01 


0.12 


0.69E-01 


c2 


0.48E-01 


0.16E-01 


0.14E-02 


-0.15E-02 


-0.13E-01 


0.32E-01 


0.22E-01 


c3 


0.39E-02 


0.14E-02 


0.73E-02 


0.20E-02 


0.19E-02 


0.14E-01 


-0.29E-02 


c4 


-0.32E-01 


-0.15E-02 


0.20E-02 


0.34E-02 


0.20E-02 


-0.39E-02 


-0.60E-02 


c5 


-0.36E-01 


-0.13E-01 


0.19E-02 


0.20E-02 


0.15E-01 


-0.15E-01 


-0.24E-01 


c6 


0.12 


0.32E-01 


0.14E-01 


-0.39E-02 


-0.15E-01 


0.18 


0.12 


c7 


0.69E-01 


0.22E-01 


-0.29E-02 


-0.60E-02 


-0.24E-01 


0.12 


0.15 



TABLE X. Covariance matrix for the seven scaling function fit parameters after applying the SciBooNE production measure- 
ments in the F-S fit. 





p8.89 


0k 


& K. prod 




(GeV/c) 


(rad) 


(mb) 


E v = 0.35 GeV 


1.52 


0.213 


9.05 ± 0.62 (6.9%) 


E v = 0.65 GeV 


2.07 


0.127 


10.32 ± 0.62 (6.0%) 


E v = 0.90 GeV 


2.45 


0.103 


9.87 ± 0.58 (5.9%) 


Kaon Sweet Spot 


2.80 


0.106 


8.37 ± 0.49 (5.9%) 


HE Events 


4.30 


0.055 


4.57 ± 0.27 (5.9%) 



TABLE XL Differential cross section values for various kine- 
matic points as in Table [V] but including in the F-S fit the 
SciBooNE production measurement for neutrino and antineu- 
trino. 



general result. Nevertheless, the SciBooNE K + neutrino 
rate measurement can be directly applied to MiniBooNE 
analysis as a constraint on the electron and muon neutri- 
nos from K + decay. Electron neutrinos from K + decays 
are one of the important background in the ty, to v e oscil- 
lation search. Understanding this background will result 
in a reduction of the systematic uncertainty in the Mini- 
BooNE oscillation analysis. 



Scaling Fits 


1.2 < P K 89 (GeV) < 5.5 






Value 


Error 




cl 


11.37 


0.93 




c2 


0.87 


0.13 




c3 


4.75 


0.09 




c4 


1.51 


0.06 




c5 


2.21 


0.12 




c6 


2.17 


0.43 




c7 


1.51 


0.40 


Input Error 


Alcshin 


1.11 


0.07 


0.10 


Allaby 


1.07 


0.06 


0.15 


Dekkers 


0.87 


0.06 


0.20 


Vorontsov 


0.54 


0.04 


5.00 


Abbott 


0.78 


0.07 


0.15 


Eichten 


1.03 


0.06 


0.15 


xVd.o.f. (no /) 


2.28 


(d.o.f. = 119) 





TABLE XIII. Results for the F-S fits as in Figure llXl but for 
the F-S fit results including the SciBooNE rate measurement, 
d.o.f. indicates here degree of freedom and "no /" means no 
correction factor applied. 



^-mode 


v-mode 


Combined 
v+v mode 


K+ Rate 0.94±0.05±0.11 


0.54±0.09±0.30 


0.88±0.04±0.10 



TABLE XII. K + rate measurement results relative to the MC 
beam prediction for the neutrino, antineutrino, and combined 
neutrino and anti-neutrino samples. Errors include statistical 
and systematic errors. 



This SciBooNE K + rate measurement has been in- 
cluded in a version of the F-S fit and the best fit results 
for the parameters including the normalization for the 
data sets is reported in Table IXIIII The covariance ma- 
trix is reported in Table IXIVI and correlation matrix is 
displayed in Figure [T3] Table IXVI lists the differential 
cross section values for kaon production at several kine- 
matic points. 

In order to apply the SciBooNE constraint to the Mini- 




FIG. 13. Correlation between the fit parameters as in Fig- 
ure [10] but for the F-S fit results including the SciBooNE rate 
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cl 


c2 


c3 


c4 


c5 


c6 


c7 


cl 


0.84 


0.47E-01 


0.39E-02 


-0.31E-01 


-0.36E-01 


0.12 


0.69E-01 


c2 


0.47E-01 


0.16E-01 


0.14E-02 


-0.14E-02 


-0.13E-01 


0.32E-01 


0.22E-01 


c3 


0.40E-02 


0.14E-02 


0.73E-02 


0.20E-02 


0.19E-02 


0.14E-01 


-0.33E-02 


c4 


-0.31E-01 


-0.14E-02 


0.20E-02 


0.34E-02 


0.20E-02 


-0.38E-02 


-0.61E-02 


c5 


-0.36E-01 


-0.13E-01 


0.19E-02 


0.20E-02 


0.15E-01 


-0.15E-01 


-0.24E-01 


c6 


0.12 


0.32E-01 


0.14E-01 


-0.38E-02 


-0.15E-01 


0.18 


0.12 


c7 


0.69E-01 


0.22E-01 


-0.33E-02 


-0.61E-02 


-0.24E-01 


0.12 


0.16 



TABLE XIV. Covariance matrix as in Table IXl but for the F-S fit results including the SciBooNE rate measurement. 





p8.89 

(GeV/c) 


0k 
(rad) 


prod 

(mb) 


E v = 0.35 GeV 


1.52 


0.213 


9.12 ± 0.62 (6.8%) 


E v = 0.65 GeV 


2.07 


0.127 


10.39 ± 0.62 (6.0%) 


E u = 0.90 GeV 


2.45 


0.103 


9.94 ± 0.58 (5.8%) 


Kaon Sweet Spot 


2.80 


0.106 


8.43 ± 0.49 (5.8%) 


HE Events 


4.30 


0.055 


4.60 ± 0.27 (5.8%) 



TABLE XV. Differential cross section values as in Table IXII 
but for the F-S fit results including the SciBooNE rate mea- 
surement. 



BooNE neutrino event prediction, one needs to consider 
the K + kinematic regions that contribute to the two sam- 
ples. 



Figure [S] shows that the kinematic region of K + mesons 
that produce background v e events in MiniBooNE and 
Figure [8] shows the regions that contributes to the Sci- 
BooNE rate measurement. While there is a large over- 
lap between the SciBooNE and MiniBooNE regions, the 
MiniBooNE region extends to somewhat lower K + mo- 
menta. Using MC studies combined with the covariance 
matrix associated with F-S fit, we have quantified the 
increased uncertainty associated with extrapolating the 
SciBooNE measurement to the lower MiniBooNE region 
and found that the error on the constrained electron neu- 
trino interaction rate should be increased by a factor of 
1.5. This increases the uncertainty for the MiniBooNE 
electron neutrino event rate prediction from the mea- 
sured SciBooNE uncertainty of 12% (as reported in Ta- 
ble IXIII) to a total error of 18%. (The associated covari- 
ance matrix given in Table IXIVI should also have all of 
the elements multiplied by (1-5) 2 = 2.25) After applying 
the new SciBooNE constraint, the MiniBooNE predic- 
tion for electron neutrinos from K + decay is reduced by 
only 3% but the uncertainty is reduced significantly by a 
factor of three from previous estimates because both the 
rate and cross section uncertainty is reduced [23j . 



VII. SUMMARY AND CONCLUSIONS 

The F-S parameterization given in Eq.[5]has a theoret- 
ically motivated form that takes into account low beam 
momentum production thresholds from exclusive chan- 
nels in contrast to many other models. For example, the 
S-W parameterization does not have the proper scaling 
properties or expected behavior for the Xf < regions. 
Also, extrapolations using data at much higher beam mo- 
mentum appear to have difficulty describing lower mo- 
mentum K + production measurements. 

The F-S parameterization describes the K + produc- 
tion data well for beam momentum in the range of 8.89 
to 24 GeV/c. Fits involving different experimental data 
sets have been performed and show good agreement with 
the experimental data as shown in Figure [5] where the 
data have been scaled by the normalization factors given 
in Table HID The normalization values (except for the 
Vorontsov data) arc in good agreement within the 10% 
to 20% uncertainties quoted by the experiments. 

The F-S fits including the full covariance matrix can be 
used to predict K + production for low beam momentum 
neutrino experiments such as the BNB at 8.89 GeV/c. 
The overall uncertainty from the fit is about 7% and is 
consistent with the combination of the experiments with 
~15% uncertainties. The fits also give the dependence on 
produced K + kinematics in angle and momentum, which 
is important for accurate neutrino flux predications using 
magnetic horn focusing devices. 

A cross check of the F-S parameterization using neu- 
trino data from the SciBooNE collaboration measure- 
ment reported in Ref. 22j confirms the accuracy of the 
model at low primary beam momenta and its validity as 
a better representation of K + production with respect to 
the S-W model. The F-S parameterization derived from 
the low energy kaon production experiments including 
this SciBooNE production constraint should therefore be 
a good representation of K + production for low energy 
neutrino beam simulations. We, therefore, suggest that 
the parameters shown in Table HXl be used along with the 
covariance given in Table IXl 

We wish to acknowledge the MiniBooNE and Sci- 
BooNE Collaboration for the use of their neutrino simu- 
lation programs and the National Science Foundation for 
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